Overview

We begin with deterministic compartmental models, which are usually straightforward to specify. With the right software, they are straightforward to solve. These models find application in virtually every scientific field. In public health, we see them most often in infectious disease epidemiology, cancer biology, pharmaco-kinetics, and healthcare operations research.

We will use examples from Brauer, Castillo-Chavez, and Castillo-Chavez (2001), Vynnycky and White (2010), Otto and Day (2011), Institute for Disease Modeling (2019).

The focus here is on specifying systems of ordinary differential equations, not on solving those systems analytically.

Overview: ODEs made easy

Write the components (“compartments” or “types” or “states”) in the time-dependent system you want to model. Examples: \(S(t)\), \(I(t)\), \(R(t)\)
Write the time-derivatives of the components as a function of time, and other compartments. Example: \(dS/dt\).
Write the initial conditions for all the components. Example: \(S(0)=99\), \(I(0)=1\).
Simulate the dynamics of the system

ODEs are most interesting when there are interactions between compartments.

Examples of deterministic models

The most boring model (part 1)

Consider the most boring model for a quantity \(y(x)\) \[ \frac{dy}{dx} = \beta \] with initial condition \(y(0)=\alpha\). The rate of change of \(y\) as a function of \(x\) is constant.

What function of \(t\) has slope \(\beta\) for every \(x\)? The linear function \[ y(x) = \alpha + \beta x . \] Have you seen this before?

The most boring model (part 2)

Recall the linear “regression” model \[ y = \alpha + \beta x + \epsilon \] where \(\epsilon\) is a mean-zero error term.

If \(x\) is a treatment, and we regard \(y(x)\) as the potential outcome when treatment is set to \(x\), then the “treatment effect” is the change in \(y\) induced by a unit change in \(x\), \[ \frac{dy}{dx} = \beta \] This is one reason that some people (economists) think of causal effects as derivatives of outcome functionals.

Birth process

Let \(y(t)\) be the number of bacteria in a dish at time \(t\). Suppose at time \(t=0\), there are \(y(0) = y_0\) bacteria. Each bacterium is immortal, and divides into two with rate \(\lambda>0\). The number of bacteria at time \(t\) is governed by the differential equation \[ \frac{dy}{dt} = \lambda y(t) \] with initial condition \(y(0)=y_0\). The solution is \[ y(t) = y_0 e^{\lambda t} \] You can verify this by computing \(dy/dt\). The dynamics of this system obey “exponential growth”. The speed of this growth is governed by \(\lambda\).

Birth process

Death process

The “death” process is the same, but we let \(\lambda=-\mu<0\). Then \[ \frac{dy}{dt} = -\mu y(t) \] with initial condition \(y(0)=y_0\), which has solution \[ y(t) = y_0 e^{-\mu t} \] and \(y(t)\) decreases/decays “exponentially” toward zero.

Logistic growth I

Recall the exponential growth model \[ \frac{dy}{dt} = \lambda y(t) \] with initial condition \(y(0)=y_0\). Its solution is \(y(t) = y_0 e^{\lambda t}\). As \(t\to\infty\), \(y(t)\to\infty\). This unbounded growth may be unrealistic. For example, resource constraints may prevent cells from dividing indefinitely.

Suppose \(K\) is the “carrying capacity” of the medium, and \(\lambda\) is the intrinsic gowth rate as above. Consider the model \[ \frac{dy}{dt} = \lambda y(t)\left( 1 - \frac{y(t)}{K}\right) \] with initial condition \(y(0) = y_0\). This is just like the pure birth model, except there is a negative quadratic term.

Logistic growth II

\[ \frac{dy}{dt} = \lambda y(t)\left( 1 - \frac{y(t)}{K}\right) \]

\(y(t)\) is constant when \(dy/dt=0\), which occurs whe either \(y(t)=0\) or \(y(t)=K\). Therefore

If \(y(t) < K\), then \(dy/dt>0\) and the population size increases toward \(K\).
If \(y(t) > K\) then \(dy/dt<0\) and the population size decreases toward \(K\).

In fact, the solution can be written down precisely: \[ y(t) = \frac{K y_0}{y_0 + (K-y_0)e^{-\lambda t}} \]

Logistic growth

Logistic growth: SI model

Consider a population of individuals where every individual is susceptible or infected with an infectious disease. Let \(S(t)\) be the proportion of susceptibles, and let \(I(t)\) be the proportion of infectives. Then the system is characterized by \[ \frac{dS}{dt} = -\beta S(t) I(t) \] \[ \frac{dI}{dt} = \beta S(t) I(t) \] where \(S(t) + I(t) = 1\). We can rewrite this model as a logistic growth model with \(K=1\) and \(\lambda=\beta\) and \(y_0=I(0)\) above.

Note that because the system only goes in one direction, the first and second equations are mirror-images.

SIR model

Consider a population of \(N\) individuals where every individual is susceptible or infected or recovered. Let \(S(t)\) be the number of susceptibles, \(I(t)\) be the number of infectives, \(R(t)\) the number of recovered. Then the system is characterized by

\[ \begin{aligned} \frac{dS}{dt} &= -\beta S(t) I(t) \\ \frac{dI}{dt} &= \beta S(t) I(t) - \gamma I(t) \\ \frac{dR}{dt} &= \gamma I(t) \end{aligned} \]

with initial condition \(S_0\), \(I(0)\), \(R(0)\) and conservation equation \(S(t)+I(t)+R(t)=N\).

SEIR model

Now suppose there is a latent period following infection that delays onset of infectiousness. Individuals in this latent state are exposed, but not yet infectious.

\[ \begin{aligned} \frac{dS}{dt} &= -\beta S(t) I(t) \\ \frac{dE}{dt} &= \beta S(t) I(t) - \delta E(t) \\ \frac{dI}{dt} &= \delta E(t) - \gamma I(t) \\ \frac{dR}{dt} &= \gamma I(t) \end{aligned} \]

SEIR model with vital dynamics

Let’s incorporate births and deaths and model the absolute number of individuals in each compartment instead of the proportion. Let \(N(t) = S(t) + E(t) + I(t) + R(t)\). Let \(\mu\) be the per-person birth rate, and let \(\nu\) be the per-person death rate.

\[ \begin{aligned} \frac{dS}{dt} &= \mu N(t) - \nu S(t) -\beta S(t) I(t) \\ \frac{dE}{dt} &= \beta S(t) I(t) - (\nu + \delta) E(t) \\ \frac{dI}{dt} &= \delta E(t) - (\nu + \gamma) I(t) \\ \frac{dR}{dt} &= \gamma I(t) - \nu R(t) \end{aligned} \]

SEIRS (loss of immunity)

\[ \begin{aligned} \frac{dS}{dt} &= -\beta S(t) I(t) + \sigma R(t) \\ \frac{dE}{dt} &= \beta S(t) I(t) - \delta E(t) \\ \frac{dI}{dt} &= \delta E(t) - \gamma I(t) \\ \frac{dR}{dt} &= \gamma I(t) - \sigma R(t) \end{aligned} \]

Example: SIR models for policy

Krebs et al. (2019)

Example: HIV care pipeline

\[ \begin{aligned} \frac{dD}{dt} &= \lambda_D I(t) - (\lambda_L + \mu) D(t) \\ \frac{dL}{dt} &= \lambda_L D(t) - (\lambda_A + \mu) L(t) \\ \frac{dA}{dt} &= \lambda_A L(t) - (\lambda_S + \mu) A(t) \\ \frac{dS}{dt} &= \lambda_S A(t) - \mu S(t) \\ \end{aligned} \]

In one-way queueing models you can think of the time spent in each compartment as the inverse of the rate of exit.

Adapted from Gonsalves et al. (2017).

Example: HIV care pipeline

Gonsalves et al. (2017)

Pharmacokinetic models

Let \(G(t)\) be the amount of drug in the GI tract, \(B(t)\) in the bloodstream, and \(O(t)\) in the organ of interest. Let \(I(t)\) be the known input dynamics of drug administration (oral, IV, etc).

\[ \begin{aligned} \frac{dG}{dt} &= I(t) - \lambda_G G(t) \\ \frac{dB}{dt} &= \lambda_G G(t) - \mu_B(t) - \lambda_T B(t) \\ \frac{dO}{dt} &= \lambda_T B(t) - \mu_O O(t) \end{aligned} \]

The \(\lambda\) parameters represent rates of transition to the next compartment, and \(\mu\) parameters represent rates of clearance.

Incorporating heterogeneity: SIR with two susceptible types

You can often incorporate heterogeneity by splitting compartments.

\[ \begin{aligned} \frac{dS_1}{dt} &= -\beta_1 S_1(t) I(t) \\ \frac{dS_2}{dt} &= -\beta_2 S_1(t) I(t) \\ \frac{dI}{dt} &= (\beta_1 S_1(t) + \beta_2 S_2(t)) I(t) - \gamma I(t) \\ \frac{dR}{dt} &= \gamma I(t) \end{aligned} \]

Solving ODEs numerically

General-purpose ODE solver: `ode`

Lots of ODEs can’t be solved analytically. For (well-behaved) systems that have a unique solution, it is often possible to solve for \(y(t)\) numerically, without having any calculus-style analytic insight into the structure of the solution. There are very fast computer algorithms for doing this.

We will focus on methods implemented in standard software packages.

library(deSolve)

you will need to install this.

For help, look at ?deSolve

Example: SIR solution

library(deSolve)
sir <- function(time, state, parameters) {
    with(as.list(c(state, parameters)), {
        dS <- -beta * S * I
        dI <- beta * S * I - gamma * I
        dR <- gamma * I
        return(list(c(dS, dI, dR)))
    })
}

Example: SIR solution

init <- c(S = 1 - 1e-06, I = 1e-06, R = 0)
parameters <- c(beta = 1.4, gamma = 0.14)
times <- seq(0, 100, by = 0.1)
out <- as.data.frame(ode(y = init, times = times, 
    func = sir, parms = parameters))
head(out)

##   time         S            I            R
## 1  0.0 0.9999990 1.000000e-06 0.000000e+00
## 2  0.1 0.9999989 1.134469e-06 1.494107e-08
## 3  0.2 0.9999987 1.287017e-06 3.189078e-08
## 4  0.3 0.9999985 1.460074e-06 5.111940e-08
## 5  0.4 0.9999983 1.656401e-06 7.293359e-08
## 6  0.5 0.9999980 1.879127e-06 9.768099e-08

Example: SIR solution

plot(times, out$I, ylab = "Population proportion", 
    xlab = "Time", type = "l", bty = "n", 
    ylim = c(0, 1), col = "red", lwd = 2)
lines(times, out$S, col = "green", lwd = 2)
lines(times, out$R, col = "blue", lwd = 2)
legend(50, 0.5, c("S(t)", "I(t)", "R(t)"), 
    lty = 1, col = c("green", "red", "blue"))

References

Brauer, Fred, Carlos Castillo-Chavez, and Carlos Castillo-Chavez. 2001. Mathematical Models in Population Biology and Epidemiology. Vol. 40. Springer.

Gonsalves, Gregg S, A David Paltiel, Paul D Cleary, Michael J Gill, Mari M Kitahata, Peter F Rebeiro, Michael J Silverberg, et al. 2017. “A Flow-Based Model of the HIV Care Continuum in the United States.” Journal of Acquired Immune Deficiency Syndromes 75 (5): 548–53.

Institute for Disease Modeling. 2019. “SEIR and SEIRS Models.” https://institutefordiseasemodeling.github.io/Documentation/general/model-seir.html.

Krebs, Emanuel, Benjamin Enns, Linwei Wang, Xiao Zang, Dimitra Panagiotoglou, Carlos Del Rio, Julia Dombrowski, et al. 2019. “Developing a Dynamic Hiv Transmission Model for 6 Us Cities: An Evidence Synthesis.” PloS One 14 (5). Public Library of Science: e0217559.

Otto, Sarah P, and Troy Day. 2011. A Biologist’s Guide to Mathematical Modeling in Ecology and Evolution. Princeton University Press.

Vynnycky, Emilia, and Richard White. 2010. An Introduction to Infectious Disease Modelling. Oxford University Press.

Deterministic models

Forrest W. Crawford

Overview

Overview: ODEs made easy

Examples of deterministic models

The most boring model (part 1)

The most boring model (part 2)

Birth process

Birth process

Death process

Logistic growth I

Logistic growth II

Logistic growth

Logistic growth: SI model

SIR model

SEIR model

SEIR model with vital dynamics

SEIRS (loss of immunity)

Example: SIR models for policy

Example: HIV care pipeline

Example: HIV care pipeline

Pharmacokinetic models

Incorporating heterogeneity: SIR with two susceptible types

Solving ODEs numerically

General-purpose ODE solver: `ode`

Example: SIR solution

Example: SIR solution

Example: SIR solution

References

Deterministic models

Forrest W. Crawford

Overview

Overview: ODEs made easy

Examples of deterministic models

The most boring model (part 1)

The most boring model (part 2)

Birth process

Birth process

Death process

Logistic growth I

Logistic growth II

Logistic growth

Logistic growth: SI model

SIR model

SEIR model

SEIR model with vital dynamics

SEIRS (loss of immunity)

Example: SIR models for policy

Example: HIV care pipeline

Example: HIV care pipeline

Pharmacokinetic models

Incorporating heterogeneity: SIR with two susceptible types

Solving ODEs numerically

General-purpose ODE solver: ode

Example: SIR solution

Example: SIR solution

Example: SIR solution

References

General-purpose ODE solver: `ode`